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ABSTRACT 


This thesis concerns developing a program that takes an 
aerial photograph, and a set of Digital Terrain Elevation 
Data CDTED) that is defined over the area of the photograph, 
and generates a synthesized view that represents what a 
Camera would see from a different location. The elevation 
data points are grouped into triangular panels that are 
projected to the reference image by three dimensional 
transformation equations. Shading for the synthesized image 
is determined from the reference image. The pixels of the 
reference image that Fall within a triangular panel are 
collected and averaged. When a new observer location is 
Selected, the panels are projected to the new synthesized 
image plane. A z2-buffer approach and a polygon Fill 
algorithm were used to remove hidden surfaces of the 
synthesized view. 

This program is tested on both artificial and real data. 
Other characteristics and performance measurements of the 
program are also analyzed here. The quality of the syn- 
thesized image from real data was affected by the low 
resolution of the terrain elevation data, and yielded less 
desirable results than could be expected of a higher 


resolution terrain model. 
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he INTRODUCTION 


A. COMPUTER IMAGE GENERATION FROM AERIAL PHOTOGRAPHY 

The main objective of this study was to develop a 
program that takes a digital photographic image and a file 
of terrain elevation points defined over that image as 
input, then produces as an output a synthesized perspective 
View. The synthesized view 1s a rotated 3-dimensional C3D) 
perspective representation of the original photographic 
image. The main application of this study is to generate a 
different perspective of a terrain model. This may be used 
to generate different views that a pilot of an aircraft 
could expect Following different Flight paths through the 
same area. Further study may make it feasible to generate 
synthesized images Fast enough to simulate a real time image 
display of a Flight for a mission briefing or to be used as 
a training aid. Another application could be for training 
men on optically guided missiles. With high resolution 
images, a Flight path through a battlefield could be 
simulated that would have all the visual characteristics of 
an actual flight without the expenditure of a live missile. 
The generation of a shaded image as a 3D picture provides 
unique problems for 3D graphic displays. The data which 
comprises a photographic image consists Of am aanaguen 


pixels, each of which has a defined grey level or shade. 
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There are @S6 different levels of grey that may be assigned 
Po. Cet neS ScuGl Cameerms taking a photographic 
perspective image and a ec-dimensional C20) array of eleva- 
tions defined ina grid covering the area of the image as 
inputs. The grid of elevations, called the terrain model, 
is geometrically related to the photographic image through a 
perspective projection transformation that equates the world 
coordinates of the elevation points to the object coor- 
dinates of the image. A synthesized image from a different 
observer location is then generated. The new synthesized 
view should approximate what would be seen by a camera from 
the new observer position. 

The differences between the original and synthesized 
images will be affected by the resolution of both the 
photographic image and the terrain models. Higher resolu- 
tion of the original models will result ina closer approx- 
imation in the synthesized image. Another complication or 
ambiguity arises when details which should show up in the 
sunthesized view were not present in the original image. A 
method must be devised so that it can Fill in areas which 
become visible in the synthesized image that were hidden in 
the original reference image. Ihe solution to this Ridden 
surface problem is further addressed in the discussion of 


the grey scale referencing algorithm.CRef. 1] 
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B. TERRAIN ELEVATION AND PHOTOGRAPH AS INPUTS 

In this study a high aititude aerial image of Moffett 
Field, California was used as the original reference photco- 
graph. The photographic image was supplied by the Defense 
Mapping Agency CDMA) and had a resolution of approximately 1 
meter per pixel. The terrain model corresponding to the 
reference image was provided as Digital Terrain Elevation 
Data CDOTED) by the DMA and consisted of elevation points 
taken every second of a degree change in latitude and 
longitude. This gives an approximate resolution of 30 
meters per elevation point ina north-south direction and 23 
meters per elevation point in an east-west direction. 

The synthesized view was restricted to a northerly 
direction which simulates an aircraft Flying from south to 
north with the image plane perpendicular to the direction of 
Flight. fo allow for different FIlIQGNt petteciswvan te 
require development of an algorithm that would provide for 
rotation of the image coordinates which is beyond the scope 
oF this study. The main idea is to generate a synthesized 
view that is rotated From the original photograph by 
approximately 93O° and explore the concepts of the algorithms 
required to do this. Although speed was not a major issue, 
the size of the terrain model was limited to a SO x SO grid 
arcay, or @S00 data points, to minimize the time For 


synthesized image generation. 


ie 


fees lodr l init.| SSUES 
ieee Ore Scale “Referencing 

To determine the grey scale values of the pixels 
that make up Our synthesized view, the terrain model is 
First divided into triangular panels. The vertices of the 
triangular panels are then mapped into the original DMA 
image using the perspective projection transformation that 
projects georectangular coordinates into reference image 
coordinates. The pixel grey scale values that fall within 
the projected triangular panel are then averaged. The 
average grey scale value is permanently assigned to that 
particular panel. When the synthesized image is constructed 
the triangular panel is mapped to the new image view and 
Filled with the assigned average grey scale value. In this 
way the sample of pixels that Fall within the triangular 
panel are mapped from the original reference image to the 
synthesized image. This method of mapping the triangular 
panels to the synthesized image also solves the problem of 
assigning a grey scale value to hidden surfaces of the 
reference image because they are automatically assigned the 
value of the surrounding pixels. Since the average grey 
scale value for a triangular panel is dependant upon the 
resolution of the terrain model, the grey scale value 
assigned to the hidden surface will also be affected 


eel. 14. 


if 8. 


The smaller the triangular panels are, the smaller 
the area that must be collected and averaged in the Original 
image. This means a much better synthesized view can be 
constructed that contains more of the attributes of the 
Original image. For this reason the resolution of the 
terrain and reference image model is very critical to 
obtaining an accurate synthesized view. Using the resolu- 
tion of the terrain and image models used in this study, the 
approximate number of pixels that must be averaged in the 
reference image for each triangular plane would be 1/2 (30m 
x e€3miC1 pixel/m) = 345 pixels. This is very coarse and 
does not allow for optimal generation of the synthesized 
view. 

©. Hidden Surface Elimination 

There are surfaces that may be discernible in the 
reference image but become hidden in the synthesized view. 
The z-buffer algorithm was used to accomplish the hidden 
surface elimination. The z-buffer is an array that 
contains the depth or distance to the observer location for 
each pixel that is to be visible in the synthesized tegen 
As each triangular plane is mapped to the synthesized view 
the location and depth of each pixel within the plane is 
determined. The depth of the pixel to be written at a 
certain location is compared to the depth of any pixel seas 
may have been previously written to the same location. le 


the depth or distance of the new pixel to the observation 
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Daiwa Ses NGGten tinal ENG ppevious pixel, its depth is 
written to the z-buffer and the grey scale value of the 
Bt wel s Splacedminto, the synthesized image. If the depth is 
MaGcder Tome eeabing OCecuEs and a neu’ pixel is obtained in 
the process. The z-buffer works very closely with the next 
mi@nichnm to be considered. (CRef. €, pp. e65-267) 
3. Polugon Fill Algorithm 

Screen coordinates are generated for the three 
vertices of each triangular panel as 1t is mapped to the 
synthesized image. Screen coordinates are designated as IA 
and JA values with the IA values representing the columns 
and the JA values the rows. The location, IA,JACO,O), 
designates the upper left hand corner of the screen and the 
maximum screen coordinate, IA,JACS12,51e2), the lower right 
hand corner. An active edge list CAEL) is generated by 
computing a line between each of the translated vertices. 
For each line, the IA coordinate corresponding to the 
maximum JA value of the line, the amount IA changes for each 
one unit step of JA, and the total span of JA are stored 
into the AEL. The three lines generated for each translated 
triangular plane will form another closed triangle. By 
using the parameters stored in the AEL the location of 
pixels enclosed by the translated triangular plane can be 
determined, and the corresponding array points within a 
Frame buffer are changed from a Otoa il. 


Moet. 2, po. 76-79) 


hee: 


After all of the enclosed pixels have been marked 
Within the frame buffer, the buffer is scanned row by row. 
If a value of 1 is Found, then the depth is calculated for 
that point and compared with the depth value stored in the 
Z-BUrT er. If the depth value is smaller, that pixel is 
located closer to the observer location and the greu Seana 
value for that pixel is written to the synthesized image 
File. As can be seen the Fill and hidden surface algorithms 
work together to generate the new image. The implementation 
of these algorithms are explained in further detail later. 

In Chapter II there will be a discussion of basic 
photographic geometry to develop an understanding of the 
transformation equations necessary to map object coordinates 
into image coordinates and for image plane rotation. 
Chapter III will detail program considerations based on 
image and elevation data as well as the algorithms used to 
generate the synthesized view. Chapter IV will discuss 
possible ways to improve the transformation program and 
discussion of topics for possible further study. An outline 
of the program is contained in Appendix A that gives a short 
discussion of each subroutine as to itS purpose, input and 
Output, modules that called, and modules that reference the 
Subroutine. Appendix B contains the entire 3D transforma- 


LLeneot Salam. 
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PieeeeenemeGRAPHIC GEOMETRY 


A. BACKGROUND 
To understand many of the concepts used in this study, a 
basic background in photographic geometry is presented. The 
relationship between the image space and object space is the 
basis For many of the equations that nelp to generate the 
reference and synthesized images for visual display. The 
Objective of this chapter is to present the concepts that 
allow the transformation of 30D objects to a 2D image and the 
parameters evolved. 
Tl. Perspective and Parallel Projection 

Fe parallel projection is a projection in which the 
projection lines from the object to the image plane never 
converge. When an object is viewed by parallel projection, 
its size would never change as the camera is moved closer or 
Further away. In contrast, a perspective projection has all 
the projection lines from an object converge to a perspec— 
eave Center. A perspective projection imitates how we see 
things. An example would be a picture of railroad tracks. 
The tracks would appear to become closer together when 
Further away from the observation point. In a parallel 
projection the tracks would be the same distance apart along 
Pema tilwne lengua. ERGE. 3, po. 1323-134) 

Since a camera 1S generally designed to photograph a 


rather large area, it involves perspective projection. [nhs 


i> 
Ni 


Camera view represents what an observer would see standing 
at the same location, and the images generated are perspec- 
tive images. This means that the equations used to trans- 
Form the object space into the image space must be perspec- 
tive transformations. 

ae Image Coordinate Space 

The image plane is the plane of the photograph to 

which the object points are mapped. It has a 2D coordinates 
System to which each point of a 3D object is translated to 
CFig. @.1). The indicated principal paint) (1 Pe ese 
center of the image plane and has the coordinates of 
Cx,y,0). The x, y, and z axis represent a right handed 
plane and the perspective center CL) lies along a line 
parallel to the z-axis; then a perpendicular line is drawn 
From L to the image plane. The point at which this perpen- 
dicular line intersects the image plane is called the 
principal point Co). This offset of the printipal sage 
From the IPP is compensated for in the transformation 
equations by xo and yo. The focal length of the plane is 
defined as the distance from the principal point to the 
perspective center. For the image plane coordinate system, 
each object point CA) is graphed to a corresponding image 
plane point Ca) located at (Cxa,ya,0O), and the perspeceive 
center or focal point is located at (xoa,yo,f). EReEh oe 


bjoe] 364 
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In generating a synthesized view, the perspective 
center may be placed at any location desired with reference 
to the image plane. It is therefore desirable to select a 
Pome aloamng cme 2Zeaexis such that xo and yo become ©. This 
Will decrease the number of calculations required in 


Generating the synthesized view. 


a Object Coordinate Space 


The observer location and each object point position 
is described in world coordinates, called georectangular 
coordinates, of X, Y, and 2. The center of the earth is 
gevem aS (0,0,0), the Z-axis points directly to true north, 
the X-axis points to the intersection of O° Latitude and O° 
Longitude, and the Y-axis the intersection of O° Latitude 
Seems ees LemGgitude. [he Sprincipal or Focal point that 
wOUld describe the observer location in georectangular 
coordinates is XL, YL, and @L. Each object point CA) 
located in the object space is identified by XA, YA, and ZA 


Semeswioume ii Figure @€.ce. ERef. 3, p. 136) 


B. IMAGE PLANE ROTATION 

To align the x, y, 2 coordinates of the image plane to 
the desired viewing direction requires rotation about the %, 
Y, and @Z@ axis of the georectangular coordinate system. ie, 
general to transform one 3D coordinate system requires a 
Poeagi~ MULE plication of the Form A= CHMIB. The A repre- 
sents a vector in the image space with x, y, and z2 coor- 


Ginates, and B is a vector in the georectangular coordinate 


Ih 
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~L-perspective centerc 





Fig. ¢.1 Image Plane Coordimates (Ref. 305573 o 





Fig. 2.2 Object Plane Coordinates (Ref 2) peo 
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Heeeneoneiee Asean 2 COmpements. This may be written as 


Ce amt eee 


x sly, ol) ie Xx 
y = qd ef Y Ce. 
Z Ga. Z 
where 
a Bio Ce Cosxx CosxyY Cosx2 
QS de Ff - CosyX CosyY CosyZ Ce ie.) 
eo ia es CoszxX CoszyY Cos22z 
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This maps the vector xX, Y, @ to the image space x, y, Z. 
The M matrix is derived from the definition of the direction 


OF a vector A given by 


A A A A 
pee fesse i + Cese j + Cosr k 
| 
oe, Ay RAS, 
ea To I pace RNG Meera: 4 C253) 
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A 
where a ie and k are unit vectors of the particular 


coordinate system in which vector A is contained. The 
Quantities «, 8, +r are the angles that vector A makes with 
the x, y, and z-axis respectively. Since there are three 
vector components in the X, Y, 2 system, each one must make 
its own transformation into x, y, and z and thereby forming 


meeeae x aS Matrix GF M. Ic translate x, y, z into xX, Y, 2 


eal 


the inverse of the ™ matrix is taken giving B = Chl A. 
CRef. 3302. 2S 

If we can define three orthogonal vectors in georect— 
angular coordinates as R, S, and T that would describe the 
desired viewing position of our image plane, the 1 matuaeaee 


easily derived by 


Ss Su a2 


ee 
Rex Ry RZ 
Sy ——, a a Cant) 


FRY IR] IR] 
ae Ty Dez 
Fe 
Generally the image plane rotation 1s expressed in omega 
Cw), phi c0), and Kappa €k). The w is the rotation about 
the X-axis, the @ is rotation about the Y-axis, and k is 
abOut. the Z-axis: If we rotated First about the X-axis By 
an angle of w radians, the terms of the M matrix would 
equate as Follows 
Cosxx ™— Coas(07*).= 1 
CosyY = CosCw) 
CosuZd = Cos(90* =v essen cw 
CoszyY,*Cestu + > S025 = ssi nc 
Cosz2Zz = Coslw) 


All other terms equate to cos 90° = O, therefore the fM 


matrix becomes 
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ote 


1 0 0 
em 0 ROS Cleon el Carraro 


O -—-Sintw) CosCw) 


Similarly a rotation @ about the Y-axis produces 


Cost @3 0 Saji Cel, 


i 0 sh 0 Cee. 


Sincd) 0 Costa) 


and for a rotation k about the Z-axis 


CosCk) Sinmck J 0 
NM = —-Sintk) CoscCk) 6 Co) 
0 0 1 


By multiplying all three matrices together we derive the 


SveGall 1 transform inw, 0, and k, 


Seas Guvbasecs ees tnISinck) + SimCwISinctICosck) 
1 = Bak (Cues seat Cle) SasGnveost ki y= Simcw)SintdsSinck) 
Se wiew.) = SintwICost od) 


Seu ostmem)  — CoscwISincdJe€osl ks 
Simon e2osck J e+ CoestwISintc da) Sinlk) 
CosCwiCosc( a) 
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This is the general form of the matrix that maps the 


georectangular coordinates to the image plane. CRef. 3, pp. 


|-7 -600) 
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The general form of the transformation matrix is used to 
initially map the elevation terrain model into the original 
reference image. The w, @, and kK were supplied with the 
Original DMA photograph and represents the rotation of the 
reference image plane with respect to the georectangular 
coordinates at the time the picture was taken. When we 
generate the synthesized view the image plane must be 
rotated to the desired viewing angle of the observer 
Cnortherly direction in particularj. fhis means that &@ new 
w, 9, and k must be calculated, or by defining new image 
plane coordinate axes in terms of the georectangular 
coordinates, we can calculate the terms of the M matrix 
directly using the preceding equations. I[his is discussed 
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In this chapter an in-depth analysis of the original 
image and terrain data is presented. [his includes how the 
data was referenced, the size of the data Files, how the 
data was used, and how the elevation and image data compared 
With one another. The process of translating from the 
object space coordinates to image space coordinates and then 
to screen coordinates 1s considered, and the equations used 
are given. 

The referencing, Fill, and 2-buffFer algorithms also are 
discussed in detail. How the data generated from these 
algorithms is used and put together to produce the syn- 
thesized view will be presented. Any problems that were 
encountered and the eventual solutions will be discussed in 


the appropriate section to which they pertain. 


A. REFERENCE IMAGE DATA 

The picture of Moffett Field supplied by the Defense 
Mapping Agency CDMA), was 4999 by 997 pixels in size and 
came with both a left and right image. Only the left image 
was used to generate the perspective views in this study. 
Fach individual pixel within the original image is desig- 
mated by coordinates I and J, where |! is the pixel column 


and J is the scan row. The geographic northeast corner has 


eo 


the I, J coordinates of €0,0), and the southWest Corner 
CEsls|7/ pess/Sise) 

The image display devices used were capable of dis- 
playing images that were only Sle By Ste pixetse goose 
therefore, the original image was divided into appropriate 
blocks suitable for viewing called frames. Each frame 
contains an image that is Sle by Sle pixels. The coor- 
dinates of each Frame has a Four digit I Frame and 
J Frame value, and is further identified as a leftwor Gea 
CL or RI image. The [I Frame and J Frame coordinates are 
designated in multiples of Sle which define the column and 
row location of each Frame. There are 10 frame columns and 
10 Frame rows with assigned coordinates of OOOO through 
4608. To identify a particular frame of the DMA image one 
First designates whether it is a Left or Right image and 
then give the |! Frame and J Frame coordinates. AS an 
example LOS1le10cec% would designate a Left image from the 
second column and third row. The first Frame LOOOOOQOO 
Starts in the southeast corner and the last Frame LYSb084608 
1s in the northwest corner. The disparity between the 
starting location of the frame coordinates and the I and J 
coordinates of the original image must be compensated for in 
the equations that are used to determine the location of 


individual pixels within a frame image as shown later. 
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1. Object to Reference Image Iransformation 


VeopemebD eGeaPol thers SCONVertEed From its 3D xX, Y, 
and 2 georectangular coordinates into the eD x and y image 
coordinates using the A = CMNIB equation discussed earlier. 
The vector from the perspective center to each object point 
me Getimed buy CXA=XLD, CYA=YL), and CZA-ZL). This vector is 
mapped into the reference image plane coordinates of x, y, 
and z. Since every vector in the image plane is directed 
From the perspective center or Focal point to the image 
point Ca), the z2 coordinate value is constant and equal to 
the negative of the focal point ¢€-f) of the camera. Using 
these parameters the equation becomes 
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where K is a scale factor. From this transformation the 


Following equations are obtained 
aoe ee oe GAN ete CYASYLI + ¢ CZA-¢2LIJ (3.2) 
Peto as Meo CXA-k es) + 2 CYA-YLI SF (CZA-ZL)9 €3.3) 
oe Meme Xk) tain CyYA-yY¥L) + 1 CZA=ZL9I ¢€3.4) 
Dividing the last equation into the first two and rear- 


ranging yields 
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where x and y are the e€D image plane coordinates fret ee 
6 fr i ed cag a 
The original parameters of the M matrix were 
calculated from the w, @, and k that were given by the DMA 
With the original photograph. These represent the physical 
position of the image plane in relation to the georect- 
angular coordinate system at the time the picture was taken. 
The Focal point was also supplied and is particular to the 
camera that was used to take the original photograph. When 
the synthesized image is generated, the image plane is 
Oriented to a position for the desired viewing angle, which 
means that the parameters of the M™ matrix will change and 
must be recalculated. The steps used to determine the 
desired orientation of the image plane coordinates will be 
discussed later in this chapter. 
a. Reference Image to Screen Coordinate Iransformation 
Once the x and y image coordinates have been 


calculated they are translated to I and J values of the 


eS 


ranainiage  lniseisa decomelished using the affine 
transform which represents a @€D into @D cocrdinate transfor- 


Mationmn. The equation that accomplishes this is derived from 
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where Cl and Ce are the values that translate the image 
plane origin to the [ and J coordinate system origin. They 
are calculated by setting [I and J to O and solving for x and 
y. To get the desired transformation of the image coor- 
dinates to I and J coordinates the inverse transform is 


meanen. LRef.3, p.S933 
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Agaiumeenie Grigimas {, kK; i.°9m; El, and Ce were 
Supplied by the DMA with the original image. Equation 3.8 
represents any general 2D into 2D transformation, and it is 
used both to translate the original image plane coordinates 
into the |! and J coordinates of the reference image and to 
transform the image plane coordinates of the synthesized 
View into screen coordinates of IA, and JA. 

The screen coordinates were assigned the parameters 
IA, which represents the columns, and JA, which represents 


the rows. The point IA,JAC1,1) is mapped to the upper left 


eo 


corner of the screen and IA, JACSle, Sle) is Phe toweqeasenee 
corner. 

The screen coordinates represent the location of a 
pixel within a frame. Io convert I and J origina leeeome 
dinates to IA and JA screen coordinates one needs to know 
the particular frame one is working in and compensate for 
the difference in the starting location of the Frame 
coordinates and the original image coordinates. This is 
accomplished by using the following equations, 

IA = Ct -— Tle Prame? Care) 
JA C4999 —- J Prame sae C3222@>? 

The I Frame and J Frame values must therefore be 
Given to determine the screen coordinates within a desired 
Frame. To allow Flexibility in determining which frame 
image would be used to extract the reference image, an 
interactive input of the I_ Frame and J Frame coordinates was 
appropriate. 

3. Synthesized Image Plane Rotation 

For the synthesized views that were generated, the 
affine transform parameters Cl, Ce, j, k, 1, andm that 
would map the newly rotated image plane into the screen 
coordinate system were selected. Since the image planes in 
the synthesized views were oriented for an observer looking 
north, the image plane coordinates were generated by cal- 
culating the z-axis, which points south, fram twa georee. 


tangular coordinate vectors calculated From two different 


30 


tercain Gate points along the same longitude line. By 


taking the difference Detween the X, Y, and 2 coordinates a 


third vector was formed that defined the image plane z-axis. 
The image plane y-axis was calculated by using only one of 
the terrain data points used to calculate the z-axis. By 
taking the negative of the georectangular coordinates of the 
terrain data point, a vector is produced that points 
downward through the center of the earth. This was the image 
plane y-axis. Once the z and y axes are calculated, the 
Cross product of y cross 2 was used to calculate the x-axis. 


Figure 3.1 demonstrates the resulting image plane. 






Synthesized 
Image Plane 
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Fig. 3.1 Synthesized Image Plane Coordinates 


This image plane coordinate system, from the 
observers perspective at (€0,0,-Ff), would have the x-axis 
Gemmcing left, the y-axis pointing down, and the z-axis 


POinting directly toward the observer. The screen 


Si 


coordinates have the IA axis pointing right and the JA axis 
pointing down. Therefore, the new values of C1 ana Geemeee 
IA, and JA equal ta O were as follows: 

Cli = Maximum assigned x-image value C Si 

So aT Cale 
which aligned the image plane x,yCo,o) to the screen 
coordinates of I[A,JACO,0). The j, k, 1, and m Valuss Ghee 
transformation matrix were selected to scale the image plane 
to the screen. 

Having determined the three synthesized image plane 

coordinate vectors in terms of georectangular vectors, it is 
relatively easy to generate the M natrix parameters using 
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were S, R, and T are the geocrectangular coordinate vectors 
of the x, y, and z axes of the rotated image plane. This 
defines the new transformation matrix that will be used to 


generate the synthesized view by mapping the georectangular 


vectors of the terrain data into the new image plane. 
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B. TERRAIN DATA 


The Digital Terrain Elevation Data CDTED) supplied by 
the DMA came as a rectangular grid of elevation data points. 
Fach elevation data point was recorded as an integer value 
im meters above sea level. If a particular elevation point 
was unknown, it was assigned a value of -3e767 to assure 
that it would not be confused with any valid elevation data 
points. The rectangular terrain grid listed elevation 
points every one second of a degree change in latitude and 
longitude. The southwest corner of the terrain grid was 
defined as being located at 37° ee’ 47” N. latitude and 
melee OS O33 Sue teongittde. ) From this reference point the 
elevation data was laid out in 210 rows by 239 columns. The 
rows represented lines of constant latitude and the columns 
lines of constant longitude. With this information the 
northeast corner of the terrain grid was calculated as being 
mgeaeee at SJ7mees 17” N. latitude and =lee O1’O"%” wW. 
Hoengictude. 

1. Data Verification 

The First problem was to compare how accurately the 
elevation data matched up with the original image data. 
This required taking specific elevation points that were 
known to match specific image points, then translating the 
georectangler coordinates of those elevation points to the I! 
and J coordinates for comparison to the original image. The 


w, ©, and k for the M matrix to transform georectangular to 


| 


image plane coordinates and the parameters for the affine 
transform from image plane to original image I and J 
coordinates were supplied by the DMA with the original image 
data. 

To select specific elevation points for compar fae 
required finding a method of distinguishing unique elevation 
patterns that could match specific objects in the images 
The technique used was to visually display the elevation 
data as an image. The elevation image File was produced by 
assigning grey scale values to each elevation data point 
with the lower elevations receiving the darker shades and 
the higher elevations the lighter shades. When the cle wa. 
tion image was displayed as shown in Figure 3.e, a distinct 
highway pattern emerged from which three intersections could 
reasonably be distinguished. The three intersection points 
Cshown as a, 5, andc iin Fig. 3.e€) correspond to elevated 
roads that crossed one another in the original image. 

Once the elevation points were selected for com- 
parison, the approximate row and column of the elevation 
data corresponding to the center of the intersections was 
determined. From the reference point of the terrain grid 
there is a 1 second change in latitude and longitude for 
each row and column which allowed the calculation of the 
latitude and longitude for each of the three reference 
elevation points. The latitude and longitude for each point 


was converted to georectangular coordinates using a 


eh 


conversion program, then transformed to |! and J coordinates 
Peemeteenseovieed UNA parameters as explained earlier. 
Fach of the three elevation points mapped to within 10 
Bixels of the original image in Both the I and J coordinate 
Sirections. This equates to less than 1 second error in 
latitude and longitude which was deemed precise enough to 


establish the correlation between the terrain and image 


data. 





Fig. 2.2 Wetevation image 
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By selecting a smaller area of the DIED data, a more 
Gistinct picture could be studied. An intersection used to 
verify image and elevation correlation was selected to be 


used as the reference image. A smaller rectangular set of 
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terrain data that would map to the intersecraon wee sea 
small section of the Surrounding area, was extracted from 
the reference terrain grid. This smaller set of terrain 
data was taken From rows 71 through 793 and columns 172 
through 183 of the terrain neces. 

To verify that this set of elevation points wou 
appear like the desired reference image, a line drawing 
illustrated in Figure 3.3 was created using a commercial 
graphics program called MOVIE.BYU CRef. %t). this proegean 
Can generate a connected line drawing From a set of eleva- 
tion data and allows rotation as well as magnification of 
the drawing. To assure a sharp visual contrast, the 
elevation data was magnified by a Factor of 10 and the 
drawing rotated to a useful viewing angle. 

The results were not as desired which gave an early 
indication that the resolution of the terrain data may not 
be adequate enough to generate a synthesized view that would 
be a close approximation of the reference image. This was 
Further verified when the synthesized view was produced at a 
later time. An artificial set of image and elevation data 
was used to verify that the transformation program func- 
tioned as desired before generating synthesized views of the 
original reference image. The artificial elevation points 
mapped into the artificial image exactly way as the original 
elevation data would map to the original reference image. 


The artificial reference image, shown in Figure 3.4, 
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was Of a small square structure resting on top of a much 
larger square structure. Selecting a lagpge Cheer reewea 
artificial reference image would minimize the effects of the 
low resolution of the elevation data. This produced more 
reasonable synthesized images that demonstrated the perseeae 
tive transformation more accurately. The results of the 
transformation will be discussed further on in this chapter 
after considering some of the algorithms used to generate 


the synthesized view. 


EC, *SPECIE LE ASG tins 
ioe Image Referencing Algorithm 

Due to the resolution mismatch between the reference 
image and the terrain data, a smaller subset of the terrain 
model was selected. This would allow the transformation of 
smaller and more distinct images that could show the effects 
of the program better. The desired terrain elevation points 
are extracted from the larger reference terrain grid By 
interactively selecting the proper rows and columpsmesaec 
define those particular points. The program allows up to a 
SO by SO array of elevation points to be extracted. T[his 
Size was chosen to limit the time required tc gepeaaweec 
synthesized view. From the rows and columns, the latitude 
and longitude is tabulated For each point, and is then 
converted tc georectangular coordinates. The georectangular 
coordinates are written to a File in row order From letra 


right and from bottom to top. The First Set SF coeradina-e 
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match the southwest elevation point, and the last set of 
coordinates the northeast elevation point. 

The georectangular coordinates of each terrain data 
point is translated to !, and J original image coordinates 
using the camera parameters supplied by the DMA. They are 
further processed to [A and JA screen coordinates using the 
I Frame and J Frame of the desired reference image. The IA 
and JA coordinates derived from the elevation points will 
now correspond to the IA and JA ccordinates of the reference 
image. Any four adjacent elevation points will define a 
Pectangle which is divided into two triangular panels by 
defining a line conmmecting two opposite corners. Once the 
triangular panels of the elevation points are defined in 
terms of IA and JA coordinates, the corresponding pixel grey 
scale values of the reference image that Fall within the 
same screen coordinates of the triangular panel are col- 
lected and averaged. As seen in the example of Figure 3.5, 
the calculated average grey scale value is placed into a 
mee ealongGg with the three specific elevation points that 
make up the triangular panel. This procedure is repeated 
until the entire terrain grid has been processed. Each 
triangular panel represents a sample or average grey scale 
Value of the reference image. 

Drepmemec lace euce, longitude, and elevation of the 
new observation point is input into the program, a new XL, 


YL, and @lL is calculated that corresponds to the new 
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perspective center. AS explained carlier, Gne tnageeoeee 
is rotated to the desired viewing angle, which Changes c7oen 


matrix parameters as well. 
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Fig. 3.5 Triangular Panel Data File Steeecece 


With these new parameters for the perspective transformation 
equations, the georectangular coordinates of the terrain 
Grid are once again run through the perspective transfor- 
mation, and the new IA and JA coordinates are calculated for 
each point. These new IA and JA coordinates represent the 
transformed elevation points as seen from the new observer 
location. The next step is to take the same three elevation 
points that formed a triangular panel in the original 
transformation, map them into the synthesized image file 


using the new IA and JA coordinates, and then Fill them in 


a0 


with the assigned grey scale value. The mapping and filling 
Deecess 15 Explained in the next section. 
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In the perspective transformation of the terrain 
Reraeinea the SYneneSiZead wmage plane, some of the tri- 
angular panels become partially aor entirely hidden by other 
panels. To determine which pixels will be visible and 
therefore written to the synthesized image and which pixels 
are hidden, the z-buffer algorithm was used. 

When the new observer location is input into the 
program, the XL, YL, and 2@L georectangular coordinates are 
tabulated. Using the georectangular coordinates previously 
Calculated for each of the terrain grid elevation points, 
the distance or depth from the observer location to the 
elevation data points can be calculated by the following 
equation 

teen Sea be GaetCCxbewe + CYL-YI* +C€ZL-2Z)?*) Se oie we 

The depth of the pixels within a triangular panel 
were calculated by using the normalized plane equation that 
defines the plane of the triangular panel. The normalized 
Dlane equation in 3D space is given by 

eo + cz et — j Ge Peas. 
To solve for the coefficients a, ©, and c, the three eleva- 
tion points that specify a triangular panel are used, and 
the x, y, and z are replaced with IA, JA, and the Depth of 


each elevation point. If the three points are CIAI1, JAI, 
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Depthli), CIAe, JAe, Depthed, and C1AS, Jas) Eeptnseeee enone 


matrix form we have the following 
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IAl1 JAI Depthl a = 
IA2 JAC Depthe b ™ =i C3.169 
IA3 JAS Depth3 e ml 
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Solving For the coefficients a, 5, and c we have CRefl ae 
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The inverse of the elevation point matrix is 
determined by calculating the adjoint, which is the trens— 
pose of the cofactor matrix, and multiplying each term by 
the reciprocal of the determinant CRef. 5, p. A-15S3 79a 
depth of each pixel within a tranformed triangular panel can 
now be determined by using the IA and JA of the pixel in the 
Following equation. 

Deptipe]- Cl + a CIA -495 CJA we C2) Be 
Each triangular panel will have different a, Bb, andc 
coefficients, therefore the pixels within each triangular 
panel will have a different depth than those from another 
panel. 

The z-buffer is a two dimensional array that is the 


same size as the synthesized image of Sie by Sle pixels. 
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Tennehnelamanoe Im COSChGimates OF a pixel is determined, the 
depth is calculated and then compared to any previously 

mere ten Geochim the z-buffer at that same coordinate 
Pmeocatlon. If the depth is smaller, it means that the pixel 
is closer to the observer and would cover the previously 
Written pixel. The depth of the new pixel will then replace 
that of the previous pixel. The assigned grey scale value, 
which 1s determined from the projected triangular panel in 
which the new pixel is contained, is written to the syn- 
thesized image at the calculated [A and JA coordinate 
location. If the new pixel depth is larger then the 
previous pixel depth, that means the new pixel is farther 
away From the observer and would be covered by the previous 
pixel. In this case no action is taken and the next pixel 
is processed. This procedure is continued until all the 
pixels within each translated triangular panel has been 
Sgoacessed. (CRef. €, pp. e&S-c67) 

The IA and JA coordinates of the pixels in the 
Synthesized image are calculated from each transformed 
triangular panel using an active edge list CAEL), J Bucket, 
Seemobame butter. the IA and JA screen coordinates of the 
three elevation points that define a triangular panel are 
used to generate the parameters of the AEL. Three lines are 
Formed that connect each of the three translated elevation 


points and are identified as line 1, 2, and 3. A line is 


op a, 


defined between any two points From which you can determine 
the [A associated with the maximum JA of Che Cue spores 

TA_INCPT = ¢€ IA value of the maximum JA ) C3. Se 
How much [A changes for a one step change in JA iS givens 

IACPoint 2) = 1eG@eenn eae 
DELTA IA = eae iaediiemeeeeeniaiatan (2. 2am 
JACPoint 2) {Je Ce etneaele 

which 1s the inverse slope of the line between the two 
points. The span of JA between two points is given by 

DELTA.JA.™. JACPointe) = WACeeameae) C3 eum 

For each line, the IA coordinate that corresponds eee 

the maximum JA value of the line, the amount IA changes for 
each one unit step of JA, and the total span of JA is 
determined and stored into the AEL. After the line para- 
meters are placed into the AEL, the line identification 
number is put into the J Bucket, which iS a one dimensional 
array of size Sle, at the same location as the maximum JA of 
the line. In this way the J_Bucket acts as a pointer to the 
lines stored in the AEL. If a previous line number has 
already been written to that location, then the line ident- 
ification number of the new line is placed into the AEL and 
1s linked to that line already residing in the J_Bucket. An 
example is shown in Figure 3.6 where line #1 is referenced 
by the J_Bucket and line #e is linked to line #1. The IA 
and JA coordinates of the pixel locations to be mappeacyes 


the synthesized image are contained within the three lines 


ie 


whose parameters are contained in the AEL. If a line is 
horizontal, the IA and JA coordinates of that line are 
located between the other two lines and therefore does not 


feed CO be Weltten to the AEL. CRef. e, pp.75-78) 
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Fig. 3.6 Active Edge List Storage 


The J Bucket is scanned from itsS maximum to minimum 
coordinate value. If the J Bucket contains a line pointer, 
the parameters for that line are retrieved from the AEL as 
well as those of any line it is linked to. Since we have 
defined a triangle there will always be two lines to work 
with. From the parameters of the two lines, the IA coor- 
dinates that fall between them is calculated for each JA 
scan line. A scan line is all the columns CIA coordinates) 
along a particular row CJA coordinate). As the JA coor- 
dinate is decremented by one, From maximum to minimum, the 
J Bucket is checked to see if a new line has been added, and 
then the corresponding IA for each line is determined for 
that JA value. Those IA values and any that fall between 


them are mapped to the frame buffer. 


oS 


The frame buffer is a Sle by Sle array that is 
initialized to 0. As the IA and JA values are mapped into 
the frame buffer the coordinate locations matching enema 
and JA values are changed from O0 to 1. This process 
continues, and each time the JA scan line is decremented, 
the DELTA JA parameters for each of the two lines are also 
decremented and the J Bucket is checked. If the J Bucket 
contains another line pointer, then the line it references 
Will replace the line whose DELTA JA has decreased to O. 
When Finished, the frame buffer will have recorded the IA 
and JA coordinates for every pixel location within the 
translated triangular panel. The frame buffer is then 
scanned, and if a particular IA and JA coordinate location 
contains a value of 1, then the depth for a pixel located at 
those same coordinates is calculated and compared to 
previously tabulated depths in the z2-buffer as explained 
San Ier. 

This process is known as a polygon fill routine that 
utilizes the z2-buffer algorithm for hidden surface ellmina- 
Elan. If any part of a triangular panel, after going 
through the perspective transformation, should map outside 
the synthesized image coordinate boundary, the entire panel 
is discarded and the next panel is processed. This 1S not a 
Satisfactory solution and could have been corrected by 
implement:ng a clipping routine that would allow parscuas 


triangular panels to be mapped to the synthesized image. 
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However, due to time constraints, a clipping routine was not 


Pneen@eeaeeemmeo the 3H transformation program. 
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Using the artificial reference image from Figure 3.4 and 
an artificial terrain grid that mapped to the same image, 
two synthesized views were generated. the first synthesized 
view shown in Figure 3.7 was From an observation point 
Beca Capaumepen cco MSO "Ne laticude. -lee, O01’ SS” WwW. 
longitude, and an elevation of 110 meters. Figure 3.8 
exhibits the next synthesized view that was produced from an 
BeSesVvVation germe Of 37° 20° OY N. Latitude while main-~ 
taining the same longitude and observer height as before. 
This simulates viewing the object from further away. As 
expected of a perspective view the object appears smaller. 
The near view also demonstrates the perspective relationship 
by the apparent taper from Front to Back. A third perspec— 
tive view is depicted in Figure 3.93 from close in and at a 
Higher elevation. The observer location was from 37° 23’ 
Pee Seeeaetrcice, -lee’ Ol” SS’°’* W. Longitude, and a height 
of 160 meters. This demonstrates the visual effect of not 
displaying partial triangular panels in the image and why a 
Clipping routine is needed. 

The actual reference image is shown in Figure 3.10 from 
which the synthesized view displayed in Figure 3.11 was 
generated. In comparing the reference image to the syn- 


thesized view there is little resemblance. A closer 
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synthesized view is seen in Figure 3.12 which gives a 
partial representation of the reference image. The lack of 
detail in both shading and the depiction of features can be 
attributed to the poor resolution of the tetratmede aaa 
proximately ¢@S meter resolution) as compared to that of the 
reference image data (1 meter resolution). 

To improve the quality of the synthesized view the 
terrain data must be of a higher resolution. Having wineee 
elevation data points over a given area, the fewer number of 
reference image pixels one must collect amd average Goa 
triangular panel. The synthesized image will then maintain 
a closer approximation to the various shades of the refer- 
ence image. The physical shape of an object also suffers 
From poor terrain data resolution. The perspective trans- 
Formation forms straight lines between translated elevation 
points. This causes object distortion if the elevation 
Points do not fall exactly along the boundaries of the 
object. As an example, if a square box in the reference 
image had only one elevation point defined on its surface at 
the center of the box, then the synthesized image can not 
reproduce the corners and edges of that box, and it would 
appear distorted. For these reasons the closer the resolu- 
tion of the terrain data to the resolution of the reference 
image, the closer the synthesized view resembles the refer- 


ence image in shading and shape. 
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Pig. 3.7 Artificial Synthesized Image 1 
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Fig. 3.93 Artificial Synthesized imagem 





Fig. 3.10 Reference stage 


50 





Fig 3.11 Synthesized Image 1 
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The 3D transformation program was developed to allow 
tracing the Flow of data easily. Much of the data was 


written to Files so that it could be printed out and 


studied. This method of data storage required certain files 


to be read many times as the synthesized image was gener- 
ated. The program execution would have been Faster if the 
data had been stored in arrays that are easily passed 
between the various modules. Another consideration that 
would increase speed would be to decrease or eliminate the 
interactive input required. fhis could be aceoenpaogecea 
extracting the desired terrain data into a file and sepa- 
rating the associated reference image before program 
execution. This would require longer set up time but would 
decrease the amount of data manipulation required by the 


program and increase its overall speed. 
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wee CONCLUSIONS 
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The 3D computer image transformation From a photographic 
image was a difficult task to achieve satisfactorily. The 
main objective was to develop a program that takes a grey 
scale photographic image, a set of elevation data points 
defined over that image, and generates a rotated synthesized 
perspective view. This goal was realized, but some areas 
still need improvement. The quality of the synthesized view 
is judged on how well the grey scale values of pixels match 
those of the original image and how closely the translated 
objects resemble the desired structure. 

The resolution of the terrain model as compared to the 
resolution of the reference image data, extensively affects 
the synthesized image quality. Depending on the contents of 
the reference image, the required resolution of the terrain 
model will vary. IF the image is of open country, the 
distance between elevation points may be large and the 
result may not suffer unacceptable degradation in the syn- 
thesized view. If the image is of a city that has many 
small distinct objects such as buildings, then the resolu- 
tion of the elevation points must be higher. Given a set of 
elevation data points, the question to be resolved is how to 
make the synthesized pixels relate to the reference image 


better, and thus improve the quality of the synthesized 


Se 


Views?’ This question and alternative ways to improve speed 


and program Flexibility will be discussed. 


B. GREY SEALE Correlate 

There are a few possible methods to improve the shading 
of the synthesized view to better match that of the original 
image. The translated triangular panels form very distine— 
tive lines or boundaries Detween areas of different shades. 
It may be desirable to blend these boundaries by sampling 
the pixels along both sides, replacing those pixels with an 
averaged value. Another possibility would be tc implement a 
Gouraud or Phong shading algoritm CRef. 3, pp. Ses-aeen 
This would help smooth the shade transition across the 
boundary and result in a more smooth appearance. 

Another method to improve grey scale correlation would 
be tc develop a way to divide the triangular panels into 
smaller triangular areas before referencing the original 
image. By having smaller triangular panels, smaller areas 
of the reference image are sampled resulting ina better 
approximation in the synthesized view. 

Only the left image of a stereo photographic pair of 
images was used in this study. The right image may contain 
grey scale information of surfaces hidden in the left image 
view or vice versa. By sampling both left and right images 
and complementing them, some ambiguities may be resolved. 


These suggestions will increase the number cf calculations 


Sit 


required to generate the synthesized view but may improve 


BReMoeal tet Gf Ehe sumenesized view. 


eee Ee ROGRAM SPEED AND FLEXIBILITY 

Although speed was not a prime consideration in this 
Peeauicceuould De desirable tO generate synthesized views 
as quickly as possible. To determine which subroutines 
consumed the most CPU time the program was monitored while 
cunning the artificial reference data set. Table 1 shows 
the results of the CPU time used and the percentage of the 
total time for each subroutine called by the main program. 
If a subroutine calls another subroutine that time is 
included in the calling routines CPU time. The subroutines 
that required interactive input were not measured. There- 
Fore, the total CPU time cannot be measured precisely. 
However, the results obtained represent reasonable estimates 
and demonstrate where the Focus should be on improving 
overall speed. 

Some suggestions have already been discussed on how to 
improve the speed of the program, but other methods also 
exist. Preconditioning the input data to eleminate interac 
tive input and using arrays instead of Files were two 
methods already presented. The z-buffer is accessed several 
times during synthesized view generation. he wee 2s oUbf er 
could be implemented in hardware, the time required for 
processing of polygon Fill and hidden surface elimination 


routines would be improved. 
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TABLE 1 


CPU TINE CONSUMPTION (iN =SEseaiinoe 


SUBROUTINE CPU TIME PERCENTAGE (%) 
TER CROP 14.11 3.051 
READIMAGE 3.23 2.072 
TER INTRP 0.04 0.026 
REAL EL 6) (7 0.045 
TER DMS 0.693 0.443 
IM REFIJ 0.33 0.250 
IM REFAUG ily LIC 0.706 
AFFIN 0.04 0.026 
NEW IJ 0.80 0.513 
NODE DEPTH 0.63 0.404 
Brie 134.80 86.467 
Toros 155.9 1OGh@ 


a Ta a ee a na 


The frame buffer is used and accessed in two different 
areas of the program for each triangular panel translated. 
It may be possible to eliminate this buffer by processing 
each pixel as its IA and JA coordinates are determined, 
instead of storing that information into the frame buffer 
for later processing. Other polygon Fill routines such as 
seed fill algorithms or using fence registers may be faster 
than the edge fill routine used in this program CRef. e, pp. 
S@=S6). 

For program flexibility it would be desirable to be able 
to adjust the image plane of the synthesized view to any 
desired angle. The program limits the image plane to a 


northern direction. To make the synthesized image plane 
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adjustable would require developing a method for rotating 
the image Plane coordinates in terms of the georectangular 
coordinates. This would improve the 3D transformation 
program and extend its usefulness. Another program improve- 
ment would be the incorporation of a clipping routine to 
improve the appearance of partial synthesized views as 
discussed in previous sections. 

The ideas used to develop this program were contrived 
From fundamental concepts. Many areas were discussed that 
could improve or enhance the basic implementation of the 
Beegram. the results that were obtained are encouraging and 


eeulad Sasily be used as the basis For further study. 
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SUNERES2ZE0 Vem 
Ge: cae .2ng Treatermites 


pe eel ene Ele 


(D 


Tailed routines 


Dise 72. 


my 


Routine parameters 


None 


te. SUBRPEU TP HEIN e wa 
a. Funeticns performed 
Calculates the new [A, and JA secreen coordinates cf 
the extracted elevation data points using Che newvesesoeees 


POeCGation Gaba. 


Ee Galelene 
Nig 3b nee Gecrectangular coordinates of the 
new observer location. 
Pocus. Assigned focal length For the 


synthesized view. 


Al (Ae. SS lee. 
rey oe Oe re a re Tre affine transform cscceff iclereae 
TERJOM: The number of rows OF the ex cheaeeem 


elevation data points. 


() 


IENDN: 


JA: 


Called routines 


euee we xl CGR TEN: 


fae wemeer of columns of the 


extracted elevation data peints. 


Pies eaee! eemtaining Ene elevaticn 
data [A screen cocrdinates cf the 

SeQenesi Zea view. 

The array containing the elevation 
data JA screen cocrcinates cf the 


Sees coee. ban. 


Routine parameters 


fel cloe seme 


XIMA: 


PAESteraonu COoOeEFIiClTents of the 


Machisce Oematlom:. Matrix. 


Seance TOW coefFicrents of the 


SEanSe@namaelone Matti, 


Pita Tene eOSrbrelLencs of the 
@eetS: GemMaelOm Matrix. 

Pee eertSee oF Ene primeipal point 
cope ne ler. 


Image plane x coordinate. 


etl 


Dt oo Bike 


= 


YIMA: 


PSO. 


IP: 


Soe eee seen 


Image plane y coordinate. 
Denominator of the trans srmaeuee 
Mace i= 

Total number of extracted elevation 


data points. 


Counter. 


PUNE LONS Derr ertied 


Calculates the distance From each transformed 


elevaticn data point toa 


oe 


e 


‘@o 


mae 


Ae tbe 


LERMEM: 


ances 2. 


the new observer location. 


The georectangular coordinates of 
the new abserver location. 

The number of rows in the extracted 
elevation data points. 

The number of columns in the 


extracted elevation data points. 


Array of the distances from the 
transformed elevation data pages 


te the mew observer lLeGgakzanme 


Poutine parameters 


ae 


a eesumeutlens performed 
Determines tre hidden surfaces using a z-bufFfe 


eeeges chm. the depth is compared For each pixel at a 


Fy 


Sm@eseifried screen coordinate te determine if it is w 


fhe sunthsesized image. 


ee Gel eitis 

Te: Array containing the IA screen 
coardinates cf the transformed 
Eelevaticr data pcints. 

ya Array containing the JA screen 
See get pa ecesmeanwe ne “fans cr Ted 
elevation data points. 

Leer: Array of the distance From the 
weageeeetece (eve Jaller Gaca DOInts 
toa the new cbserver locaticn 

IMAGE : Peead USSG 6LS Construct ehe 
sunthesized image. 

ee ace Ese maAeMnaSaeer rows OF “tre extracted 
SlevaeltemiGdca Ceimcs’. 

Ee, ee the number cf columns oF the 


Brocade eee S&S Svactlion data eoints. 


Be 


0) 


(D 


my 


Ae Ope File of the georectanculer Geen 


dinates of the elevation data 


pe rAesS: 

CEG ea = 

[PrPSES: A File cortaining the suntheswazaa 
image. 


CalienG seems 


sib jet The 2Ze=pur Por 


Ailes The screen cocrdinatses anc desth cf 
the three elevation seints = Se 
define a triangular parel. 

a ee The distance of a pixel Withee 

rianguiarc panel to the new 


cbserver location. 


Celt. Gee. cesar 

Cal ween Gad: The cofactors cf the plane Sscvualues 
Mae ieee 

Ea. The determinate of the planes 


equation taba = 


Sie 


A COEF, 8 COEF. 


cial (@, |B/a) a ae 


REKE cs 


4 


PUT fd: 


4 


tw wT 
+ ‘ oes 


MIN: 


ee 


J MAX: 


K, and L: 


IPLANES: 


ie eceehicremes Of the silane 


ele (Ese steiae 


Grey scale valus sf a triangular 


Danel. 


“Iumerical @ssignaticn cf the thres 


Gloves lerprdoucwma= limes somal makes 
triangular panel. 

(oSntmenieele SGreen coordinate af 
the three elevation eeints. 
Maximum [IA screen cocrdinatse of 
three elevatican points. 

Minimum JA screen coordinate of 
three elevation poirts. 

Maximum JA screen ccordinate cf 


the three elevation peints. 


The frame buffer. 


CEuUnMeers : 


UWStae sem bes GF triangular panel 


constructed Fram the extracted 


elevation data peints. 


(0 
oh 


ct 
2p 
(D 


the 


S 


lo. SUBROE TT VE Iara 
a. Functions performed 
Constructs an edge list From the three transtermee 
Elsvaticon points. Thess edges Form the boundaries For a 
PelLesor bores ee The pixels that are determined to 


Fall within the trarsformed triengle are marked im Che wee 


buffer 
asl gee) obs. 

GSO) ee 2 ee | BEB Pe 

and; HEE Ec. The numerical designaticn cf the 
three elevation data points that 
make a triangular panel. 

IA: Array that contains the IA screen 
coordinates of the transfcrmed 
elevation data points. 

JA: Array that contains the JA screen 
coordinates of the transformed 
elevation data points. 

JS eS Maximum JA screen coordinate of the 
tnree elevation data points. 

Bp Minimum JA screen coordinates of the 
three elevation data pcints. 

= JU@eb 

PROMIE: The Frame buffer array. 
dim Cet sor oe oe 

ae ee 


es 


e, Called routine 


None. 


a) 


eee CoS a> Baa some hers 


Dalene Aleem e Cf “aalaine. 


DELTA | The amcunt [A changes Fer @ ans 
Seopa. La JF. 

el Tia The tstal JA span cf a line. 

ieee Array containing the parameters of 


the three lines of transformed 
triangular panel. 


XNODEL and 


AHOBE Ve: Designates the number of IA coor 
dinates between two lines for a given 
JA. 

ee Indicator usec to determine the 


Secs tte ewe tae lA eoor— 
dinates are countec. 

mous = Array containing the numerical 
designation of the threes elievaticn 
data points. 

Vee... =ne 

Hegre: foes esas Fes CeLSemiming <hre 
elevation point that contains the 
Highest JA value between twe 


Point S.. 
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Hens GH and 


ie Be 


4 
t+ 
4 
= 


ana 1S: 


RI 38/2) be 


Litandla- 


LEM “aneeceeri. 


Used with MNCDE! and NGZBaaes 
determining the highest JA value 


between two points. 


Counters. 

Array that contains the line number 
designators referencing the AEL. 
Used to determine the number of IA 
coordinates to be written to the 
Frame buffer. 


Counters for the 7 778422 
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Pe ae BS 


3D-TRANSFORMATION PROGRAM LISTING 
PROGRAM TRAN_3_D 


THIS PROGRAM TAKES AN IMAGE AND ELEVATION FILE AND 
Ceo reve he Tae ene Ph ERENCH = HACE AND ELEVATION FILE. 
FROM THESE FILES A SYNTHESIZED IMAGE IS PRODUCED. 


CHARACTER ELFILE*13, IMFILE*13 

BYTE IMAGE( 512,512) 

INTEGER IELEV2(50,50),IA( 2500), JA( 2500) 
INTEGER IROW,LROW, ICOL, LCOL, IENDN, IENDM, 
INTEGER IFRAME, JFRAME 

REAL X1,Y1,2Z1, FOCUS, DEPTH( 2500) 

REAL Al,A2,B1,B2,C1,C2, OMEGA, PHI, KAPPA 


CALL INPUT( IROW, LROW, ICOL, LCOL, ELFILE, IMFILE, 
IFRAME, JFRAME) 


OPEN( UNIT=1, FILE=ELFILE, STATUS='OLD' ) 


OPEN( UNIT=2, FILE='ZFIL. DAT' , STATUS='NEWw', 
ACCESS='DIRECT' , RECORDSIZE=128, MAXREC=512 ) 


OPEN( UNIT=3, FILE='XYZ. DAT' , STATUS='NEW', 


ACCESS=' SEQUENTIAL' , FORM='FORMATTED' ) 


OPEN( UNIT=4, FILE=IMFILE, STATUS='OLD',ACCESS='DIRECT', 


RECORDS IZE=128, MAXREC=512 ) 


OPEN( UNIT=20, FILE='NODE. DAT’ , STATUS='NEW' , 


ACCESS='DIRECT' , RECORDSIZE=128, MAXREC=512 ) 


OPEN( UNIT=21, FILE=' IMAGES. DAT' , STATUS='NEW', 


ACCESS='DIRECT' , RECORDSIZE=128, MAXREC=512) 


ete ROOF WOW, bROW, [1COL, LCOL, [ELEV2 ) 


IENDN=LROW-IROW+1 
IENDM=LCOL-ICOL+t+1 
CALL READIMAGE( IMAGE) 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


Ceyiiu 


CALL 


TER_INTRP( IENDN, IENDM, IELEV2) 

REAL_EL( IENDN, IENDM, IELEV2) 

TER_DMS( IROW, LROW, ICOL, LCOL, IENDM) 

IM_REFIJ( IA, JA, IFRAME, JFRAME, IENDM, IENDN) 

IM_REFAVG( IA, JA, IMAGE, IENDM, IENDN) 

AFFIN(A1,A2,B1,B2,C1,C2) 

OBS_LOC( X1,Y1,2Z1, FOCUS) 

Lee wiz FOCUS) All.A2,B81,B2,C1,C2, 
IENDM, IENDN, IA, JA) 

NODE_DPTH( X1,Y1,21,DEPTH, IENDM, IENDN) 

FILL( IA, JA, DEPTH, IMAGE, IENDM, IENDN) 


CLOSE( 1) 
CLOSE(2) 
CLOSE( 3) 
CLOSE( 4) 
CLOSE( 20) 
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CLOSE( 21) 
END 
CG KERR RRR KKK RRR KKK RRR KKK KKK KKK KKK 
SUBROUTINE INPUT( IROW, LROW, ICOL, LCOL, ELFILE, IMFILE, 
IFRAME, JFRAME) 


TROW : INITIAL ROW OF DESIRED AREA 
LROW : LAST ROW OF DESIRED AREA 
ICOL : INITIAL COLUMN OF DESIRED AREA 
LCOL : LAST COLUMN OF DESIRED AREA 
ELFILE : ELEVATION DATA FILE NAME 
IMFILE : IMAGE DATA FILE NAME 
IFRAME : I FRAME NUMBER 
JFRAME : J FRAME NUMBER 

CHARACTER BLEILE* 137 ib er ae 

INTEGER IROW, LROW, ICOL, LCOL, IFRAME, JFRAME 


QAQaQAQANAANANA 


WRITE(6,*)'INPUT ELEVATION AREA DESIRED’ 
WRITE(6,*)'ENTER MINIMUM ROW NUMBER : ' 
READ(5,35)IROW 
WRITE(6,*)'ENTER MAXIMUM ROW NUMBER : ' 
READ(5,35)LROW 
WRITE(6,*)'ENTER MINIMUM COLUMN NUMBER : ' 
READ( 5,35) ICOL 
WRITE(6,*)'ENTER MAXIMUM COLUMN NUMBER : ' 
READ(5,35)LCOL 

35 FORMAT( I3) 
WRITE(6,*)'INPUT THE ELEVATION DATA FILE NAME : ' 
READ(5,45)ELFILE 
WRITE(6,*)'INPUT THE IMAGE DATA FILE NAME: ' 
READ(5,45)IMFILE 

45 FORMAT(A13) 
WRITE(6,*)'INPUT THE I FRAME NUMBER OF IMAGE : ' 
READ(5,55) IFRAME 
WRITE(6,*)'INPUT THE J FRAME NUMBER OF IMAGE : ' 
READ(5,55)JFRAME 


55 FORMAT( 14) 
WRITE(6,%*)'***** WAIT APPROX. 1 MINUTE FOR SETUP*****'! 
RETURN 
END 

c 

G EK KKRKKKKRK KKK KK KKK KK KK KKK KKK KKK KKK KKK KKK KKK RRR RRR 

SUBROUTINE TER_CROP( IROW, LROW, ICOL, LCOL, IELEV2) 

c 

¢ THIS SUBROUTINE READS THE ORIGINAL TERRAIN GRID 

Cc AND CONSTRUCTS A SMALLER GRID OF THE ELEVATION 

C POINTS DESIRED. 

@ 


CHARACTER SWLAT*8, SWLON*8, DELLAT*4, DELLON*4, 
CHARACTER COLS*4, ROWS*4 
INTEGER IELEV1( 210,239), IELEV2(50,50), IROW, LROW 
INTEGER ICOL, LCOL, JV( 239) 

C THE VALUES OF JV,SWLON, SWLAT DELLON, DELLAT,COLS 
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xO 
20 


CPE OC CO Ie) 


COG a 2 @ 


10 
20 


AND ROWS ARE IGNORED. 
READ( 1,5)SWLON, SWLAT, DELLON, DELLAT, COLS, ROWS 
FORMAT( 1X,2(A8,2X),4(A4, 2X) ) 
DO 20 M=1,238 

READ( 1,10)JV(M),( IELEV1(N,M),N=1, 210) 
FORMAT( I6, 2x, 2016/( 8X, 2016) ) 
CONTINUE 
DO 40 N=IROW, LROW 

DO 30 M=ICOL,LCOL 

[ELEV2( M+1-ICOL,N+1~IROW)=IELEV1(N,M) 

CONTINUE 
CONTINUE 
RETURN 
END 


KERR KRREKKRERREKEKRKRERRKREKKRKEKRKREKRKRKKKRKRKRKRKKR KKK 


SUBROUTINE READIMAGE( IMAGE) 
THIS SUBROUTINE READS THE IMAGE DATA INTO AN ARRAY. 
BYTE IMAGE( 512,512) 


POeO Vint ).5 12 
ReaD 4yeec—-fR)( IMACE( IC, IR), IC=1,512Z) 
CONTINUE 
RETURN 
END 


REE EREKEKKEKEREKEKEKEKRRKREKERKREKEKEKERKEKRKRKEREKRKERKRKRKRKKRKRKRKRKRKRKRKREK 


SUrRoOUGmih TERSINIRE( JENDN; IENDM, LELEV2) 


THIS SUBROUTINE DETERMINES THE AVERAGE VALUE OF THE 
ELEVATION DATA AND ASSIGNS THAT VALUE TO ANY UNKNOWN 
DOIN es 


INTEGER IENDN, IENDM, IELEV2(50,50) 
DARAYNSXm, TEL, 1AVG/0, 0, 0/ 


DO 20 N=1,IENDN 
DO 10 M=1,IENDM 
IF( IELEV2(M,N). EQ. -32767) THEN 
Gore 10 
ELSE 
NEXT=NEXT+1 
IEL=IEL+IELEV2(M,N) 
END IF 
CONTINUE 
CONTINUE 
IAVG=IEL/NEXT 
CHANGE UNKNOWN ELEVATION VALUES TO THE 
CALCULATED AVERAGE. 
DO 40 N=1,IENDN 
DO 30 M=1,IENDM 


=) ak 


aaa exe) 


10 


20 


Qa 


IF( IELEV2(M,N). EQ. -32767) THEN 
IELEV2(M,N)=IAVG 
END IF 
CONTINUE 
CONTINUE 
RETURN 
END 
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SUBROUTINE REAL_EL( IENDN, IENDM, IELEV2) 
THIS SUBROUTINE CREATS A REAL FILE OF THE EREVATTG Ne 


INTEGER IELEV2(50,50),IENDN, IENDM 
REAL AELEV( 239) 


DO 20 N=1, IENDN 
K=0 

DO 10 M=1,IENDM 

K=K+1 

AELEV( K)=IELEV2(M,N) 

CONTINUE 

WRITE( 2, REC=N)( AELEV( INT), INT=1, IENDM) 
CONTINUE 

RETURN 

END 


RRR KEK RE KEK KKREKKRRERRKEKRKEKRRRRKEKKKEKEKREKKRKRKRREKRKREKREKE 


SUBROUTINE TER_DMS( IROW, LROW, ICOL, LCOL, IENDM) 


MQQaqqqaagrT }aQ1AQAQ 1 OQ 


THIS SUBROUTINE CONVERTS EACH ELEVATION POINT TOGeeee 

DMS EQUIVALENT. IT USES THE FACT THAT EACH BLEVATIG? 

POINT REPRESENTS A ONE SECOND CHANGE IN LATITUDE 

OR LONCITUDE FROM THE NEXT = roma. 

REFERENCE POINT: LAT.037 DEG.:22 MIN. : 47) S52@is) 2 oe 
LON. -122 DEG.:05 MIN. : 0S SEG a 

OUR ELEVATION POINTS STAY WITHIN THE BOUNDS ObmGey 

DEGREES NORTH AND -122 DEGREES WEST, SO THESE VALUES 

WILL BE ASSIGNED. 


IMPLICIT DOUBLE PRECISICNPGs=2) 

INTEGER IROW, LROW, ICOL, ECOL, TENDM, Fey ae 

REAL X,Y,2Z,ULATD, LATM, LATS, AELEV( 239) 

REAL LOND, LONM, LONS, HEIGHT 

PARAMETER( RLATM=22. , RLATS=47. , RLONM=5. , RLONS=3. ) 


LATD=37. 
LOND=-122z. 
DO 20 IL=IROW, LROW 
K=IL/60 
LATM=RLATM+K 
LATS=RLATS+( IL-K*60) 
IE( LATS. CE 60. 0 jitaan 
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PS al ON Oe, @ 


LATS=LATS-60. 0 
LATM=LATM+1. 0 

END IF 

Piles Geena) 

READ( 2, REC=I1)( AELEV( INT), INT=1, IENDM) 


I=0 
DO 


10 JL=ICOL, LCOL 


K=JL/60 
LONM=RLONM-K 
LONS=RLONS-( JL-K*60) 


IF( LONS. LT. 0. 0) THEN 


LONM=LONM=-1. 0 
LONS=LONS+60. 0 


END IF 
LONM==-LONM 
LONS=-LONS 
I=I+1 


HE TICHIT=AELEV( 1) 
CALL DMS2XY2Z( LATD, LATM, LATS, LOND, LONM, LONS, 


ee le ten) 


WRITE(3,99)X,Y,Z 
FORMAT( 3(1X,F17.7)) 
CONTINUE 
CONTINUE 
ENDFILE( 3) 
REWIND( 3) 
RETURN 


END 


KRREKKKKEKKEKK KKK KHER K KKK HKK KEK KKK KKK KKKEKKKKEKEKKKEKKRKKKREK 


SUBROUTINE IM_REFIJ( 1IA,JA, IFRAME, JFRAME, IENDM, IENDN } 


Pls eoUsKOULING, CONSTRUCIS THE I AND J COORDINATE 
DATA FOR EACH ELEVATION POINT AND THEN CONVERTS THEM 


LO 


SCREEN COORDINATES AND STORES THEM IN ARRAYS 


IA AND JA. 


IMPLICIT DOUBLE PRECISION (A-Z) 


REAL 


CMe ae rm Kare, %,24,X1,7%1,21,40, YO 


Rotel, © IMA Ady AZ, BL,B2Z,C1i,CzZ2,EOCUS 
INTEGER I,J, 1TENDM, IENDN, IA(2500),JA( 2500) 
INTEGER IFRAME, JFRAME, ITOT 


DATA 
DATA 
DATA 
DATA 
DATA 
DATA 
DATA 


ITOT= 


OMEGA, PHI, KAPPA/. 8341764, -. 4563699, 3. 0761254/ 
XO, Y0/0. 000002,0.0/ 
X1,Y1,Z21/-2693765.9,-4304520. 4,3859018. 3/ 
A1,A2/20. 11323959, -6. 022849824/ 

B1,B2/6. 016207940, 20. 10938801/ 
C1,C2/-34954. 59484, -22566. 71593/ 
FOCUS /0. 153197/ 


Rae LE NDM 


DOWvOm lL k—1, I TOL 
READ( 3,5,END=Z0)X,Y,2 


oh8. 


5 FORMAT( 3(1X,F17.7)) 
CALL PROJECT(X,Y,Z, OMEGA, PHI, KAPPA, XO, YO, X1,Y1,Z1, 
XIMA, YIMA, FOCUS) 
CALL XY2IJ( XIMA,YIMA,1,J,Al,A2, BIpEene ieee 
IA( IR) =I-IFRAME 
JA( IR)=4999-JFRAME-J 


10 CONTINUE 

20 REWIND( 3) 
RETURN 
END 
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SUBROUTINE DMS2XY2Z( LATD, LATIM, LATS, LOND, LONM, LONS, 
LONS, HEIGHT X72) 


THIS SUBROUTINE CONVERTS DMS DATA TO X,Y, AND Z 
GEORECTANGULAR COORDINATES. 


M0200 


IMPLICIT DOUBLE PRECISION (A-Z) 
REAL PHI,LAMDA,N,X,Y,Z,LATD, LATM, LATS 

REAL LOND, LONM, LONS, HEIGHT 
PARAMETER( PI=3. 14159265358793238) 
PARAMETER( C1=180. ,C2=60. ,C3=3600. ) 
PARAMETER( E_SQUARE=0. 006768658, A=6378206. 4) 


RADIAN=PI/Cl 

PHI=( LATD+LATM/C2+LATS/C3 ) *RADIAN 
LAMDA=( LOND+LONM/C2+LONS/C3 ) *RADIAN 
N=A/SQRT( 1-E_SQUARE*SIN( PHI)*SIN( PHI)) 
X=(N+HEIGHT ) *COS( PHI ) *COS( LAMDA) 
Y=(N+HEIGHT)*COS( PHI) *SIN( LAMDA) 
Z=(N*( 1-E_SQUARE)+HEIGHT) *SIN( PHI) 
RETURN 

END 


KEE KEEKREEKEEKKEKEKK KKK KRKRKRKEKREKKRKKKKKRKKKKRKRKRKRKRKRRKRKRKRKRKRKKRKKEEK 


SUBROUTINE PROJECT( X,Y,2Z,OMEGA, PHI , KAPPA) XO) oa 
Y1,21,XIMA, YIMA SOCuSe 


QQ 


THIS SUBROUTINE CONVERTS THE X,Y,4Z GEORECEANGUE 
COORDINATES TO XIMA AND YIMA WHICH ARE IMAGE 
PLANE COORDINATES. 


QQ0AaQQ 


IMPLICIT DOUBLE PRECISION (A-Z) 

REAL M11,M1i2,Mi3,M21,M22,M23,M3i) M32) Meay eee 
REAL XIMA,YIMA,X,Y,2Z,OMEGA, PH!) Kare 

REAL XO, YO,X41, Yi Zee ee 


M11=COS( PHI) *COS( KAPPA) 
M12=COS( OMEGA) *SIN( KAPPA) + 
| SIN( OMEGA) *SIN( PHI) *COS( KAPPA) 
M13=SIN( OMEGA) *SIN( KAPPA) - 

COS( OMEGA) *SIN( PHI) *COS( KAPPA) 
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OOF OOF =O? 


M21=-COS( PHI )*SIN( KAPPA) 
M22=COS( OMEGA) *COS( KAPPA) - 

SIN( OMEGA) *SIN( PHI) *SIN( KAPPA) 
M23=SIN( OMEGA) *COS( KAPPA) + 

COS( OMEGA) *SIN( PHI) *SIN( KAPPA) 
M31=SIN( PHI) 
M32=-SIN( OMEGA) *COS( PHI) 
M33=COS( OMEGA) *COS( PHI) 


DENOM=M31*( X-X1)+M32*( Y-Y1)+M33*(Z-Z1) 
XIMA=XO-FOCUS*(M11*( X-X1)+M12*( Y-Y1)+M13*(Z-Z1))/DENOM 
YIMA=YO-FOCUS*(M21*( X-X1)+M22*( Y-Y1)+M23*(Z-Z1))/DENOM 
XIMA=XIMA* 1000000 

YIMA=YIMA*1000000 

RETURN 

END 
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SUbr Cue zis( XIMA, YIMA,1,J,Al,AZ,B1,82,C1,C2Z) 


om oeemouliiee LAKES THE IMACE POINTS AIMA,YIMA 
AND CONVERTS THEM TO I,J ORIGINAL IMAGE COORDINATES. 


(MPEG re POUBLE PRECISION (A=-Z) 
READ IMaVIMaA, Al, AZ,B1l,B2,C1,Cz,DENOM 
PEE Che ld 


DENOM=A1*B2-B1*A2 

I=( (XIMA-Cl1)*B2-( YIMA-C2)*B1) /DENOM 
=-(( XIMA-Cl) *A2-( YIMA-C2)*Al)/DENOM 
RETURN 

END 


ee ee ee ee ee ee ee ee eee eee eee ee ee eee ee ee ee ee 


SUBROUTINE IM_REFAVG( IA, JA, IMAGE, IENDM, IENDN) 


‘eee hk OU LUNE CONSTRUCTS THE FILE THAT IDENTIFIES 
THE ELEVATION POINTS THAT MAKE UP A TRIANGULAR PANEL 
AND ITS ASSOCIATED SAMPLED GREY SCALE VALUE. 


IMPLICIT DOUBLE PRECISION (A-Z) 
REAL SLOPE, YINT 

BYTE IMAGE( 512,512) 

INTEGER IA(2500),JA(2500),IY,M,N, IGREY1, IGREY2,L,L1 
INTEGER IENDM, IENDN,NODE_A,NODE_B,NODE_C,NODE_D,IR 
INTEGER ICOUNT1, ICOUNT2, ITOT1, ITOT2 


DO 90 IR=1, IENDN-1 

IL=( IR-1)*IENDM 

DO 80 N=1+IL, IENDM-1+IL 
NODE_A=N 
NODE_B=N+1 
NODE_C=N+IENDM 


Sh) 


219 


60 


70 


80 
20 


exe} 


CPOG@? 02-0) 


SLOPE=( JA( NODE_C) -JA( NODE_B) )*1. 0/ 
( IA(NODE_C)-IA( NODE_B))*1.0 

YINT=1. O*JA( NODE_B)-SLOPE*IA( NODE_B) 

ITOT1=0 

ITOT2=0 

ICOUNT1=0 

ICOUNT2=0 

DO 70 M=IA( NODE_B), IA( NODE_A) 

IY=( SLOPE*M+YINT) 

DO 50 L=JA(NODE_B),IY 
ITOT1=ITOT1+IMAGE(M,L) 
ICOUNT1=ICOUNT1+1 

CONTINUE 

IF(M. LT. IA(NODE_A) ) THEN 
DO 60 L=IY+1,JA(NODE_C) 

ITOT2=ITOT2+IMAGE(M,L) 
ICOUNT2=ICOUNT2+1 
CONTINUE 
END IF 
CONTINUE 
Li=(N=( 1R=1))*2 
IGREY1=ITOT1/ICOUNT1 
IGREY2=ITOT2/ICOUNT2 
NODE_D=NODE_C+1 
WRITE( 20, REC=L1-1)NODE_A, NODE_B, NODE_C, IGREY1 
WRITE( 20, REC=L1)NODE_B, NODE_D, NODE_C, IGREY2 
CONTINUE 
CONTINUE 
RETURN 
END 
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SUBROUTINE EXT_ORIEN(Mi1,M12,Mi3s MZ M22 
M31,M32,M33,IENDM, IENDN) 


THIS ROUTINE DETERMINES THE M MATRIX PARAMETERS 
FOR ROTATION OF THE IMAGE PLANE TO DESIRED 
LOCATION FOR VIEWING IN THE SYNTHESIZED IMAGE. 


IMPLICIT DOUBLE PRECISION (A-Z) 

REAL MAGN_X,MAGN_Y,MAGN_Z,X,Y,Z 

REAL X_CORD( 2500) , Y_CORD( 2500) , Z_CORD( 2500) 
REAL X_VECX, X_VECY, X_VECZ, Y_VECX, Y_VECY, Y_VECZ 
REAL Z. VEC, 2 AVEGY Ze ee 

INTEGER IENDM, IENDN,ITOT,IP,IR 


I TOT=IENDM* I ENDN 
DO 10 IR=1,ITOT 
READ(3,5,END=20)X,Y,Z 
X_CORD( IR)=x 

Y_CORD( IR)=Y 
Z_CORD(IR)=2Z 

FORMAT( 3(1X,F17.7)) 
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10 CONTINUE 
20 REWIND( 3) 
Y_VECX=-( X_CORD(1)) 
Y_VECY=-( Y_CORD(1)) 
Y_VECZ=<-( Z_CORD(1)) 
IP=ITOT-IENDM+1 
Z_VECK=X_CORD( 1)-X_CORD( IP) 
Z_VECY=Y_CORD( 1)-Y_CORD( IP) 
Z_VECZ=Z_CORD( 1)-Z_CORD( IP) 
e USE THE CROSS PRODUCT OF Y CROSS Z TO OBTAIN 
e THE X VECTOR. 
XK VECK=(( Y_VECY*Z_VECZ)-( Y_VECZ*Z_VECY) ) 
X_VECY=( ( Y_VECZ*Z_VECX)-( Y_VECX*Z_VECZ) ) 
X_VECZ=( ( Y_VECX*Z_VECY)-( Y_VECY*Z_VECX) ) 
MAGN_Z=SORT( ( Z_VECX**2)+( Z_VECY**2)+( Z_VECZ**2) ) 
MAGN_X=SORT( (X_VECX**2)+( X_VECY**2)+(X_VECZ**2) ) 
MAGN_Y=SORT( ( Y_VECX**2)+( Y_VECY**2)+( Y_VECZ**2) ) 
M11=X_VECX/MAGN_X 
M12=X_VECY/MAGN_X 
M13=X_VECZ/MAGN_X 
M21=Y_VECX/MAGN_Y 
M22=Y_VECY/MAGN_Y 
M23=Y_VECZ/MAGN_Y 
M31=Z_VECX/MAGN_Z 
M32=Z_VECY/MAGN_Z 
M33=Z_VECZ/MAGN_Z 


RETURN 
END 
é 
C ee ee eee eee eee ee eee eee eee eee eee ee ee ee eee ee ee ee ee ee ee eee 
SUBROUTINE AFFIN(A1,A2,B1,B2,C1,C2) 
Cc 
S THIS SUBROUTINE ASSIGNS OR CALCULATES THE 
C COEFFICIENTS TO BE UTILIZED IN THE TRANSFORM FROM 
é IMAGE COORDINATES TO SCREEN COORDINATES. 
E 
IMPLICIT DOUBLE PRECISION(A-Z) 
REAL Al,A2,B1,B2,C1,C2,XIMA_MAX, YIMA_MAX, I_MAX, J_MAX 
DATA I_MAX, J_MAX/512.0,512.0/ 
C 
XIMA_MAX=1600. 0 
YIMA_MAX=1600. 0 
C1=XIMA_MAX 
OP=6)..10 
A2=0.0 
B1=0.0 
Al=-XIMA_MAX/( I_MAX*1. 0) 
B2=YIMA_MAX/( J_MAX*1. 0) 
RETURN 
END 
C 
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THIS SUBROUTINE CALCULATES THE NEW OBSERVER Xi y ae 
LOCATION FROM DESIRED LAT. AND LONG. INPUTS AS WELL 
AS PROVIDE THE FOCAL LENGr 


IMPLICIT DOUBLE PRECIS(TONeG-2Z) 
REAL LATD, LATM, LATS, LOND, LONM, LONS, HEIGHT 
REAL XA1,Y¥i,21,X%,7 7,2, Beer 


LAT D=37 50 
WRITE(6,*)'INPUT OBSERVER LATITUDE IN-MINUTES( REAL): ' 
READ(5,5)LATM 

WRITE(6,*)' -SECONDS( REAL): ' 
READ(5,5)LATS 

LOND=-122.0 

WRITE(6,*)'INPUT OBSERVER LONGITUDE IN-MINUTES(REAL):' 
READ(5,5)LONM 

WRITE(6,*)' -SECONDS( REAL): ' 
READ(5,5)LONS 

FORMAT( FS. 1) 

WRITE(6,*)'INPUT OBSERVER HEIGHT-METERS(REAL): ' 
READ(5,10)HEIGHT 

FORMAT( F6. 1) 

FOCUS=0. 015 

CALL DMS2XYZ( LATD, LATM, LATS, LOND, LONM, LONS, HEIGHT, 
MZ) 


Y1L=Y 
Z1=Z 
RETURN 
END 
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SUBROUTINE NEW_IJ(X1,Y1,21, FOCUS, AD AZ Ba sZ, Cee 
IENDM, IENDN, IA, JA) 


THIS SUBROUTINE COMPUTES THE NEW IA AND JA SCREEN 
COORDINATES FROM THE GIVEN OBSERVER LOCATION. 


IMPLICIT DOUBLE PRECISION (A-Z) 
REAL X1,Y1,Z1,FOCUS,M11,M12,M13,M21,M22,M23,M31 
REAL X0,Y0O,X,Y,Z,XIMA, YIMA,Al,A2,B1,B2,Cl,C2 

REAL M32,M33,DENOM 

INTEGER IENDM, IENDN, ITOT, IR, IA( 2500), JA(2500),I,J 
DATA. KO) YO/0. 0.08.0 


ITOT=IENDN* IENDM 
DO 10 IR=1,2500 
IA( IR)=0 
JA( IR)=0 
CONTINUE 


CALL EXT_ORIEN(M11,M12,M13,M21,M22,M23) Msi Mia2Zy oes 


IENDM, IENDN) 
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DO 20 IR=1,ITOT 
READ( 3,15, END=30)X,Y,Z 

FORMAT( 3( 1X,F17.7)) 

DENOM=M31*( X-X1)+M32*( Y-Y1)+M33*(Z-Z1) 

X IMA=XO-FOCUS*(M11*( X-X1)+M12*( Y-Y1)+M13*( 2-21) )/ 
DENCM 
YIMA=YO-FOCUS*(M21*( X-X1)+M22*( Y-Y1)+M23*(Z-Z1))/ 
DENOM 

XIMA=XIMA* 1000000. 0 

YIMA=YIMA*1000000. 0 

GALL XY213( XIMA, YIMA,I,J,A1,A2,B1,B2,C1,C2) 


IA( IR)=I 

SAC IR)=J 
CONTINUE 
REWIND( 3) 
RETURN 
END 
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SUE NOUit! Pe NOUR mer id Alp rd) Zi) DEPTH, TLENDM, [ENDN) 


THIS SUBROUTINE CALCULATES THE DISTANCE OF THE 
TRANSLATED ELEVATION POINTS TO THE NEW OBSERVER 
LOCATION. 


IMPLICIT DOUBLE PRECISION( A-2Z) 
REAL X1,Y1,21,DEPTH( 2500) 
INTEGER IENDM, IENDN,IR,ITOT 


ITOT=IENDM* IENDN 
Por loulR=),,1TOT 

READ(3,5,END=20)X,Y,Z 

FORMAT( 3( 1X, F17.7)) 

DEPTH( IR)=SORT( ( ( X1-X)**2)+((Y1-Y)**2)+((Z1-Z)**2)) 
CONTINUE 

REWIND( 3) 

RETURN 

END 
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SUBROUTINE FILL( IA,JA,DEPTH, IMAGE, IENDM, IENDN) 


Tiles SUBROUTINE DETERMINES THE HIDDEN SURFACES AND 
CONS PRUCTS THE STRANSLALTED IMAGE. A Z_BUEEER IS USED 
eae eee OMrUmED OE PERHS FROM THE OBSERVER TO THE 
GEOGRAPHIC POSITION. A PLANE EQUATION IS CONSTRUCTED 
Booleans seve POINTS. THIS BOUATION IS USED TO 
Pee eeNe THESDEPETH OF ALL POINTS WITHIN THE PLANE. 


IMPLICIT DOUBLE PRECISION( A-Z) 

BYTE IMAGE(512,512) 

REAL DEPTH( 2500) ,Z_BUFF(512,512),X1,X2,X3,Y1,Y2,Y3 
Mawes 2 DEPTH, Cli1,Ci2,C13,C21,C22,C23,C31 


Sle 


20 a7 


REAL C32,C33,DET,A_COEF , BEC@CEE Greer 

INTEGER IGREY, NODE_A, NODE_B, NODE_@) 1 INF ia 
INTEGER J_MAX, FRAME( 512,512), 1A( 2500), JA( 2500), IENDM 
INTEGER IENDN, IR, 17o7K) bees 


FIRST DETERMINE THE NUMBER OF PLANES AND INITIALIZE 
THE 2 BUPFER AND IMAG Tee. 


IPLANES=( ( IENDM-1)*2)*( IENDN-1) 
DO 20 K=1,512 
DO 10 L=1,512 
IMAGE(L,K)=0 
Z_BUFF(L,K)=0 
CONTINUE 
CONTINUE 


DETERMINE THE COEFFICIENTS OF THE PLANE EQUATTe) 
FROM THE DGRREE BLEVATION FoOthis. 


DO 70 IR=1,IPLANES 
READ( 20, REC=IR)NODE_A, NODE_B, NODE_C, IGREY 
X1=IA(NODE_A) 

Y1=JA( NODE_A) 
Z1=DEPTH( NODE_A) 
X2=IA( NODE_B) 
Y2=JA( NODE_B) 
Z2=DEPTH( NODE_B) 
X3=IA(NODE_C) 
Y3=JA( NODE_C) 
Z23=DEPTH( NODE_C) 

DETERMINE THE COFACTOR ELEMENTS 
Cl1=( (Y2*Z23)-(Y3*2Z2)) 
C12=-( (X2*Z3)-(X3*2Z2)) 
C13=( (X2*Y3)-(X3*Y2) ) 
C21=-( (Y1*Z3)-(Y3*Z1)) 
C22=(( lee) es 
C23=-((X1*Y3)-(X3*Y1)) 
C31=( (Y1*22)-( Y2*21)) 
C32=-(( X1*Z2)-(X2*21)) 
C33=( (X1*Y2)-(X2*Y1) ) 

CALCULATE THE DETERMINANT 


DET=( X1*C11)+( Y1*C12)+(21*C13) 


THE COEFFICIENTS ARE DETERMINED FROM MULTIPLYING 
THE reciprocal of the determenant with the 
transpose of the cofactors called the adjoint. 


A_COEF=-((C11+C21+C31)/DET) 
B_COEF=-( (C12+C22+C32)/DET) 

C_COEF=-( (C13+C23+C33)/DET) 

IF( C_COEE, EO, 010) GozomZO 

DETERMINE THE MAXIMUM AND MINIMUM VALUES TO TEST 


100 
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I_MAX=MAX( IA( NODE_A), IA( NODE_B), IA( NODE_C) ) 
I_MIN=MIN( IA( NODE_A), IA(NODE_B), IA( NODE_C)) 
J_MAX=MAX( JA( NODE_A) , JA( NODE_B), JA( NODE_C) ) 
J_MIN=MIN( JA( NODE_A) , JA( NODE_B) , JA( NODE_C) ) 
IF( I_MIN. LT. 1. OR. ILMAX. GT. 512)GOTO 70 
Peeve OR I IMAX, GT.512)GOTO 70 


CLEAR THE REFERENCE FRAME BUFFER AND CALL THE 
PRAME Even SUBROUTINE. 


DO 40 L=I_MIN, I_MAX 
DO 30 K=J_MIN, J_MAX 
FRAME( L,K)=0 
CONTINUE 
CONTINUE 
CALL FRAME_FIL( NODE_A, NODE_B,NODE_C, FRAME, IA, JA, 
J_MAX, J_MIN) 
DO 60 J=J_MIN, J_MAX 
DO 50 I=I_MIN, I_MAX 
LF( FRAME(I,J).EQ.1)THEN 
Z_DPTH=-( 1+( A_COEF*I)+(B_COEF*J) ) /C_COEF 
IF( Z_BUFF(1I,J).EQ. 0. 0.OR. Z_DPTH. LT. Z2_BUFE(1,J))THEN 
Z_BUFF(1,J)=2_DPTH | 
IMAGE( 1, J)=IGREY 
END IF 
END IF 
CONTINUE 
CONTINUE 
CONTINUE 
Bence d-1) 512 
PemtnOe) REC=J)CIMAGE( 1,J),1=1,512) 
CONTINUE 
RETURN 
END 
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SUBROUTINE FRAME_FIL( NODE_A, NODE_B, NODE_C, FRAME, 
IA, JA, J_MAX, J_MIN) 


THIS SUBROUTINE CONSTRUCTS AN EDGE LIST FROM THE 
THREE NODES PASSED IN THE ROUTINE. THESE EDGES 
ARE USED IN A POLYGON FILL ROUTINE USEING A FRAME 
BUFFER AND Y_BUCKET 

IMPLICIT DOUBLE PRECISION(A-Z) 

REAL I_INCPT, DELTA_I, DELTA_J, AEL( 3, 4) 

REAL XNODE1,XNODE2,DX 

INTEGER NODE_A,NODE_B,NODE_C, FRAME( 512,512), IA( 2500) 

INTEGER JA( 2500) ,J_MAX, J_MIN, NODE( 4),1S,NODE1,NODE2 

INTEGER N_HIGH, N_LOW, IT, IU, J_BUCKET(512),IR,I1,12 

INTEGER ICNT,LCNT 


agent 


10 


Zo 
30 


40 


DOwlO Ts—ay se 
J_BUCKET( IS)=0 
CONTINUE 
DO 30 Tesi 
DO 20 IR=1,4 
AEL(IS,IR)=0.0 
CONTINUE 
CONTINUE 
NODE(1)=NODE_A 
NODE( 2)=NODE_B 
NODE(3)=NODE_C 
NODE( 4)=NODE_A 
Do 46 sie—1 
NODE1=NODE( IS) 
NODE2=NODE( IS+1) 


IF( JA( NODE1). GE. JA( NODE2 ) ) THEN 


N_HIGH=NODE1 
N_LOW=NODE2 

ELSE 

N_HIGH=NODE2 
N_LOW=NODE1 

END IF 

I_INCPT=IA( N_HIGH) 


DELTA_J=( JA( N_HIGH)-JA(N_LOW) ) 


LE DELTA .2O. 0. 0) THEN 


GOTO 40 
RCSD 


DELTA_I=-( IA( N_HIGH) -IA( N_LOW) ) /DELTA_J 


END IF 
IT=JA( N_HIGH) 
IU=J_BUCKET( IT) 
IF( IU. EQ. 0) THEN 
J_BUCKET(IT)=IS 
ELSE 
AEL( IU,4)=IS 
END IF 
AEL( 1S;1)=12INGEr 
AEL(IS,2)=DELTA_I 
AEL(IS,3)=DELTA_J 
CONTINUE 
IT=J_BUCKET( J_MAX) 
IU=INT( AEL( IT, 4)) 
XNODE1=AEL( IT, 1) 
XNODE2=AEL(IU,1) 
DX=XNODE1-XNODE2 
IF(DX. LE. 0.0) THEN 
I1=NINT( XNODE1) 
I2=NINT( XNODE2 ) 
ELSE 
I1=NINT( XNODE2) 
I2=NINT( XNODE1) 
END IF 
DO 50 Ik=1 2 


LOZ 


FRAME( IR, J_MAX)=1 
50 CONTINUE 
ASC 1 ts )=AnL( 1G oe) 20 
Sew eiU os y=Anl( [Ups )=1, 0 
ICNT=J_MAX-1 
LCNT=J_MIN 
Pom Oels ren, LOND, —1 
IF( J_BUCKET( IS). EQ. 0) THEN 
XNODE1=XNODE1+AEL( IT, 2) 
XNODE2=XNODE2+AEL( IU,2) 
ELSE IF(AEL(IT,3). LE. 0. 0)THEN 
IT=J_BUCKET( IS) 
XNODE1=AEL( IT, 1) 
XNODE2=XNODE2+AEL( IU,2) 
ELSE 
[U=J_BUCKET( IS) 
XNODE1=XNODE1+AEL( IT, 2) 
XNODE2=AEL( IU, 1) 
END IF 
DX=XNODE1-XNODE2 
IF( DX. LE. 0. 0) THEN 
I1=NINT( XNODE1) 
I2=NINT( XNODE2 ) 
ELSE 
I 1=NINT( XNODE2) 
I2=NINT( XNODE1) 
END IF 
DO 60 IR=I1,12 
FRAME( IR, IS)=1 
60 CONTINUE 
Ae Mies) Ann 13 )-1. 0 
Pena Snie( 1 ,3)-1. 0 
70 CONTINUE 
RETURN 
END 
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